home *** CD-ROM | disk | FTP | other *** search
/ Graphics Plus / Graphics Plus.iso / formats / gif / dither.txt < prev    next >
Encoding:
Internet Message Format  |  1988-03-21  |  32.6 KB
  1. From: garyo@masscomp.UUCP (Gary Oberbrunner)
  2. Subject: Re: Advanced Dither Needed
  3. Date: 7 Jan 88 03:03:58 GMT
  4.  
  5. In article <3703@ames.arpa> watson@ames.UUCP (John S. Watson) writes:
  6. >    I'm looking for either references to, or code for a really 
  7. >nice color dithering algorithm.  The algorithm produces 8-bit dithered 
  8. >images that are almost indistinguishable from the original 24-bit ima
  9.  
  10. Floyd & Steinberg's dithering algorithm is described in the following paper:
  11.  
  12. Floyd, R.W. and Steinberg, L. "An Adaptive Algorithm for Spatial Gray Scale."
  13. SID 75m Int. Symp. Dig. Tech. Papers (1975), p. 36.
  14.  
  15. I don't actually have this paper, just a reference to it in Heckbert's "Color
  16. Image Quantization for Frame Buffer Display (Siggraph proceedings
  17. p. 297).  You probably want to check this paper out too.
  18. But the way F/S dither works is like this (it's so simple...)
  19.  
  20.  - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  21. First you come up with a color map - use uniform segmentation for spee
  22. Heckbert's method for quality.  Then you work from the top left to the b
  23. right of the image (pixel order), as follows:
  24.  
  25. Find the closest match in the color map to the true pixel color.
  26. Use that color map entry to represent that pixel.
  27.  
  28. Compute the r,g,b error resulting from the above approximation.
  29.  
  30. Add 3/8 of that error to the pixel to the right, and 3/8 to the pixel
  31. below the current pixel.  The remaining 1/4 goes to the pixel diagonally
  32. below.  Don't forget to clip the rgb, as well as check your boundary
  33. conditions.
  34.  
  35. And that's it - compute the next pixel!  No dither maps, no muss no fuss no
  36. bother.  And it's real fast since you can do 
  37. the 3/8 and 1/4 by shifting.
  38. On some images it can help to add some random noise; for instance if you
  39. have a smooth (flat-shaded) surface that's just below one of the color
  40. the map, the error will build up slowly until you get a line of brighter
  41. color somewhere in the surface.  A bit o' noise gets rid of that line 
  42.  - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  43.  
  44. Of course choosing the color map is the fun part. :-) :-)
  45.  
  46.                     Gary Oberbrunner
  47.  
  48. Here is a hack C program which implements both Heckbert's algorithm and F/S
  49. dither (no claims of any kind are made for this program - in fact t
  50. ended above with my name, and this is all line nois)@(#*&...
  51. No seriously, this code needs real help - it was written for a one-shot deal,
  52. and it is very inefficient in many places.  But it does work ok.
  53. The worst part is where every possible color mapping is computed - SLOW!
  54. And Remap() ain't a speed-demon either, but that's where Floyd et. al.
  55. come in...
  56.             Enjoy (watch out for possible signature at end...)
  57.  
  58. ----------------------------------------------------------------------------
  59. Remember,               -Truth is not beauty;
  60. Information is not knowledge; /    Beauty is not love;      Gary Oberbrunner
  61. Knowledge is not wisdom;     /    Love is not music;      ...!masscomp!garyo
  62. Wisdom is not truth;    ----/    Music is the best. - FZ
  63. ----------------------------------------------------------------------------
  64. ---------------------- c u t    h e r e --------------------------
  65. #! /bin/sh
  66. # This is a shell archive, meaning:
  67. # 1. Remove everything above the #! /bin/sh line.
  68. # 2. Save the resulting text in a file.
  69. # 3. Execute the file with /bin/sh (not csh) to create:
  70. #    quantize.c    (most everything in here)
  71. #    split.c        (the color-space-splitting algorithm)
  72. #    quantize.h    (structs & externs)
  73. #    std.h        (my "gotta have it" stuff)
  74. # This archive created: Wed Jan  6 21:58:18 1988
  75. export PATH; PATH=/bin:/usr/bin:$PATH
  76. if test -f 'quantize.c'
  77. then
  78.     echo shar: "will not over-write existing file 'quantize.c'"
  79. else
  80. cat << \SHAR_EOF > 'quantize.c'
  81. /*----------------------------------------------------------------------
  82.  * Color image quantizer, from Paul Heckbert's paper in
  83.  * Computer Graphics, vol.16 #3, July 1982 (Siggraph proceedings),
  84.  * pp. 297-304.
  85.  * By Gary Oberbrunner, copyright c. 1988.
  86.  *----------------------------------------------------------------------
  87.  */
  88.  
  89. #include <stdio.h>
  90. #include <fcntl.h>
  91. #include <std.h>
  92. #include "quantize.h"
  93.  
  94. static char SCCSId[] = "SCCS ID: %W% %G% %U";
  95. int NCOLS = 0;
  96. char *progname;
  97. int debug = 1;
  98. int dither = 0;                /* This turns on dithering */
  99. int fs_dither = 1;            /* This turns on Floyd-Steinberg dithering */
  100.  
  101. #define DITH 8
  102. int Dither[DITH][DITH] = {
  103. { 0, 32,  8, 40,  2, 34, 10, 42},
  104. {48, 16, 56, 24, 50, 18, 58, 26},
  105. {12, 44,  4, 36, 14, 46,  6, 38},
  106. {60, 28, 52, 20, 62, 30, 54, 22},
  107. { 3, 35, 11, 43,  1, 33,  9, 41},
  108. {51, 19, 59, 27, 49, 17, 57, 25},
  109. {15, 47,  7, 39, 13, 45,  5, 37},
  110. {63, 31, 55, 23, 61, 29, 53, 21}
  111. };
  112.  
  113. /* This is only temporarily global: */
  114.     out_pix_t  *quant_image;        /* the quantized image - cmap indices */
  115.  
  116. main (argc, argv)
  117. int argc;
  118. char **argv;
  119. {
  120.     char  *filename, *outfilename;
  121.     int    rows = 0, NCOLS = 0, cmaplen = 0;
  122.     rgb_pix_t *image;        /* rgb image storage */
  123.     cmap_t *cmap;        /* color map; rgb byte triples */
  124.     unsigned int   *hist;    /* histogram of color values: indices are
  125.                    formed from 5 bits each of r,g,b */
  126.     unsigned int   *maptable;    /* quantized color map table: indices are
  127.                    colors, values are color map indices */
  128.     int i;    /* temp! */
  129.     progname = argv[0];
  130.     if (argc != 6) 
  131.     {
  132.     dbug1("Usage: %s infile rows cols colormaplen outfile\n", argv[0]);
  133.     dbug0("Converts 24-bit rgb image to n-bit colormapped image.\n");
  134.     exit(1);
  135.     }
  136.     
  137.     filename = argv[1];
  138.     rows = atoi(argv[2]);
  139.     NCOLS = atoi(argv[3]);
  140.     cmaplen = atoi(argv[4]);
  141.     outfilename = argv[5];
  142.  
  143.     /* dynamically allocate all the various structures */
  144.     NEW(image, rgb_pix_t, rows * NCOLS * RGB_BYTES);
  145.     NEW(quant_image, out_pix_t, rows * NCOLS);
  146.     NEW(cmap, cmap_t, 3 * cmaplen);
  147.     NEW(hist, unsigned int, HIST_SIZE);
  148.     maptable = hist;        /* reuse the same storage for the map table */
  149.     
  150.     /* Perform the algorithm */
  151.     ReadImage(filename, image, rows, NCOLS);
  152.     MakeHist(image, rows * NCOLS, hist);
  153. /*    ShowHist(hist);*/
  154.     MakeCmap(hist, cmap, cmaplen);
  155.     ComputeMapping(cmap, cmaplen, maptable);
  156.     RemapImage(image, rows, NCOLS, cmap, cmaplen, maptable, quant_image);
  157.     WriteImage(outfilename, quant_image, rows, NCOLS, cmap, cmaplen);
  158.     /* Show the result */
  159. #ifdef GRAPHICS
  160.     mgiasngp(0,0);
  161.     for (i = 0; i < cmaplen*3; i+=3)
  162.     {
  163.     printf("Before LoadMap: cmap[%4d] = %3d, %3d, %3d\n",
  164.            i/3, cmap[i], cmap[i+1], cmap[i+2]);
  165.     }
  166.     LoadMap(cmap, cmaplen);
  167.     Draw(quant_image, rows, NCOLS);
  168.     mgideagp();
  169. #endif
  170. }    
  171. /* This assumes that the image is stored in packed r,g,b byte triples in */
  172. /* row-major order. */
  173. ReadImage(name, array, rows, cols)
  174. char *name;
  175. rgb_pix_t *array;
  176. int rows, cols;
  177. {
  178.     int file;
  179.     int xfercount;
  180.     file = open(name, O_RDONLY);
  181.     if (file == -1)
  182.     {
  183.     dbug2("%s: Can't open file %s for reading.\n", progname, name);
  184.     perror(name);
  185.     exit(1);
  186.     }
  187.     dbug0("Reading file...");
  188.     xfercount = read(file, array, RGB_BYTES * rows * cols);
  189.     if (xfercount == -1)
  190.     {
  191.     dbug2("%s: Can't read file %s.\n", progname, name);
  192.     perror(name);
  193.     exit(1);
  194.     }
  195.     if (xfercount != RGB_BYTES * rows * cols)
  196.     {
  197.     dbug2("File %s ended prematurely after %d bytes: continuing.\n",
  198.           name, xfercount);
  199.     }
  200.     dbug0("Done.\n");
  201. }
  202.  
  203. /* Write out the quantized image */
  204. WriteImage(name, array, rows, cols, cmap, cmaplen)
  205. char *name;
  206. out_pix_t *array;
  207. int rows, cols;
  208. cmap_t *cmap;
  209. int cmaplen;
  210. {
  211.     int file;
  212.     int xfercount;
  213.     struct FileData {
  214.     int rows, cols, cmaplen, dummy;
  215.     } fdata;
  216.     int i;
  217.     fdata.rows = rows;
  218.     fdata.cols = cols;
  219.     fdata.cmaplen = cmaplen;
  220.     fdata.dummy = 0;
  221.     file = open(name, O_WRONLY | O_TRUNC | O_CREAT, 0666);
  222.     if (file == -1)
  223.     {
  224.     dbug2("%s: Can't open file %s for writing.\n", progname, name);
  225.     perror(name);
  226.     exit(1);
  227.     }
  228.     dbug0("Writing file...");
  229.     dbug1("Header is %d bytes.\n", sizeof(fdata));
  230.     xfercount = write(file, &fdata, sizeof(fdata));
  231.     if (xfercount != sizeof(fdata))
  232.     {
  233.     dbug3("%s: Can't write file %s's header.  Wrote %d bytes.\n",
  234.           progname, name, xfercount);
  235.     if (xfercount == -1)
  236.         perror(name);
  237.     exit(1);
  238.     }
  239.     dbug1("Color map is %d bytes.\n", 3 * cmaplen);
  240.     xfercount = write(file, cmap, 3 * cmaplen);
  241.     if (xfercount != 3*cmaplen)
  242.     {
  243.     dbug3("%s: Can't write file %s's color map.  Wrote %d bytes.\n",
  244.           progname, name, xfercount);
  245.     if (xfercount == -1)
  246.         perror(name);
  247.     exit(1);
  248.     }
  249.     xfercount = write(file, array, rows * cols);
  250.     if (xfercount != rows * cols)
  251.     {
  252.     dbug3("%s: Can't write image to file %s.  Wrote %d bytes.\n",
  253.           progname, name, xfercount);
  254.     if (xfercount == -1)
  255.         perror(name);
  256.     exit(1);
  257.     }
  258.     dbug0("Done.\n");
  259.     close(file);
  260. }    
  261.  
  262. MakeHist(image, npix, hist)
  263. rgb_pix_t *image;
  264. int npix;
  265. unsigned int *hist;
  266. {
  267.     register int pix;        /* pixel index */
  268.     register rgb_pix_t *ppix;    /* pointer to the actual triple */
  269.  
  270.     /* Clear out the histogram array */
  271.     for (pix = HIST_SIZE - 1; pix >= 0; pix--) hist[pix] = 0;
  272.  
  273.     dbug0("Making histogram...\n");
  274.     for (pix = npix-1; pix >= 0; pix--)
  275.     {
  276.     if (pix % 25000 == 0) dbug1("%d \r", pix);
  277.     ppix = PIXPIX(image, pix);
  278.     /* Increment the histogram bucket corresponding to the value at pix:
  279.        Only five bits of the color value are used to increase 
  280.        clumping in the histogram.  The histogram is indexed by the five
  281.        bits each of red, green and blue all strung together. */
  282.     hist[(*ppix >> HIST_SHIFT << (HIST_QUANT_BITS * 2)) |
  283.          (*(ppix+1) >> HIST_SHIFT << HIST_QUANT_BITS) |
  284.          (*(ppix+2) >> HIST_SHIFT)]++;
  285.     }
  286.     dbug0("Done.\n");
  287. }
  288.  
  289. ShowHist(hist)
  290. unsigned int *hist;
  291. {
  292.     register int i, nmaps = 0;
  293.     printf("RED\tGREEN\tBLUE\n");
  294.     for (i = 0; i < HIST_SIZE; i++)
  295.     {
  296.     if (hist[i] != 0)
  297.     {
  298.         nmaps++;
  299.         printf("%-d\t%-d\t%-d\t= %d.\n",
  300.            i >> 10 << 3, (i >> 5 & 0x1f) << 3, (i & 0x1f) << 3,
  301.            hist[i]);
  302.     }
  303.     }
  304.     dbug1("%d different pixel values (5-bit quantized).\n", nmaps);
  305. }
  306.  
  307. /* Compute the map color for each original color */
  308. /* This should be done dynamically while remapping... */
  309. ComputeMapping(cmap, cmaplen, maptable)
  310. cmap_t *cmap;
  311. int cmaplen;
  312. unsigned int *maptable;
  313. {
  314.     register int mapcolor, d, mindist, minmap;
  315.     int orig_color;
  316.  
  317.     for (orig_color = 0; orig_color < HIST_SIZE; orig_color++)
  318.     {
  319.     mindist = BIG_DISTANCE;
  320.     for (mapcolor = 3*(cmaplen-1); mapcolor >= 0; mapcolor -= 3)
  321.     {
  322.         if ((d = DIST(((orig_color >> 10) & 0x1f) << 3,
  323.               ((orig_color >> 5) & 0x1f) << 3,
  324.               (orig_color & 0x1f) << 3,
  325.               cmap[mapcolor], cmap[mapcolor+1], cmap[mapcolor+2]))
  326.         < mindist)
  327.         {
  328.         mindist = d;
  329.         minmap = mapcolor;
  330.         }
  331.     }
  332.     maptable[orig_color] = minmap/3;
  333.     if ((orig_color & 0x3ff) == 0) dbug1("%d/31 \r", orig_color >> 10);
  334.     }
  335. }
  336.  
  337. /* Pick color map values that are close to the real ones, using the map
  338.  * computed above in ComputeMapping(). */
  339. /* Now with optional dither! */
  340. RemapImage(image, rows, cols, cmap, cmaplen, maptable, quant_image)
  341. rgb_pix_t *image;            /* original image (rgb) */
  342. int rows, cols;            /* size of the image */
  343. cmap_t *cmap;            /* color map (rgb) */
  344. int cmaplen;            /* number of entries in cmap */
  345. unsigned int *maptable;        /* mapping table from rgb5 to cmap indices */
  346. out_pix_t *quant_image;        /* quantized image */
  347. {
  348.     int npix = rows * cols;
  349.     register unsigned int pix;
  350.     register int mapentry;
  351.     int nearestcm;
  352.     dbug0("Remapping image...");
  353.     for (pix = 0; pix < npix; pix++)
  354.     {
  355.     mapentry = ((PIXRED(image, pix) >> 3 & 0x1f) << 10) |
  356.            ((PIXGREEN(image, pix) >> 3 & 0x1f) << 5) |
  357.            (PIXBLUE(image, pix) >> 3 & 0x1f);
  358.     ASSERT (mapentry >= 0 && mapentry < HIST_SIZE)
  359.     nearestcm = maptable[mapentry];
  360.     if (fs_dither)
  361.     {
  362.         register int r, g, b;
  363.         register int errr, errg, errb;
  364.         /* Compute error in current pixel */
  365.         errr = PIXRED(image, pix) - cmap[nearestcm * 3];
  366.         errg = PIXGREEN(image, pix) - cmap[nearestcm * 3 + 1];
  367.         errb = PIXBLUE(image, pix) - cmap[nearestcm * 3 + 2];
  368.         /* Propagate errors right and down */
  369.         if (pix % cols != cols - 1)
  370.         {    /* right gets 3/8 of the error */
  371.         r = (int) PIXRED(image,pix+1) +  errr * 3 / 8;
  372.         g = (int) PIXGREEN(image,pix+1) + errg * 3 / 8;
  373.         b = (int) PIXBLUE(image,pix+1) +  errb * 3 / 8;
  374.         PIXRED(image,pix+1) = clamp(r, 0, 255);
  375.         PIXGREEN(image,pix+1) = clamp(g, 0, 255);
  376.         PIXBLUE(image,pix+1) = clamp(b, 0, 255);
  377.         }
  378.         if (pix < npix - cols)
  379.         {    /* down gets another 3/8 */
  380.         r = (int) PIXRED(image,pix+cols) + errr * 3 / 8;
  381.         g = (int) PIXGREEN(image,pix+cols) + errg * 3 / 8;
  382.         b = (int) PIXBLUE(image,pix+cols) + errb * 3 / 8;
  383.         PIXRED(image,pix+cols) = clamp(r, 0, 255);
  384.         PIXGREEN(image,pix+cols) = clamp(g, 0, 255);
  385.         PIXBLUE(image,pix+cols) = clamp(b, 0, 255);
  386.         }
  387.         if (pix < npix - cols - 1)
  388.         {    /* and diagonally down gets the last 1/4. */
  389.         r = (int) PIXRED(image,pix+cols+1) + errr / 4;
  390.         g = (int) PIXGREEN(image,pix+cols+1) + errg / 4;
  391.         b = (int) PIXBLUE(image,pix+cols+1) + errb / 4;
  392.         PIXRED(image,pix+cols+1) = clamp(r, 0, 255);
  393.         PIXGREEN(image,pix+cols+1) = clamp(g, 0, 255);
  394.         PIXBLUE(image,pix+cols+1) = clamp(b, 0, 255);
  395.         }
  396.         /* Store the current pixel */
  397.         quant_image[pix] = nearestcm;
  398.     }
  399.     else if (dither)
  400.     {
  401.         register unsigned int r, g, b;
  402.         register unsigned int nr, ng, nb;
  403.         register int oppr, oppg, oppb;
  404.         int oppnrcm;
  405.         /* true rgb value for pixel */
  406.         r = PIXRED(image, pix);
  407.         g = PIXGREEN(image, pix);
  408.         b = PIXBLUE(image, pix);
  409.         /* rgb values for nearest color */
  410.         nr = cmap[nearestcm * 3];
  411.         ng = cmap[nearestcm * 3 + 1];
  412.         nb = cmap[nearestcm * 3 + 2];
  413.         /* Color as far from rgb as nrngnb but in the other direction */
  414.         oppr = min(255, max(0, 2 * (int) r - (int) nr));
  415.         oppg = min(255, max(0, 2 * (int) g - (int) ng));
  416.         oppb = min(255, max(0, 2 * (int) b - (int) nb));
  417.         /* Nearest match for opposite color: */
  418.         oppnrcm = maptable[((oppr >> 3 & 0x1f) << 10) |
  419.                 ((oppg >> 3 & 0x1f) << 5) |
  420.                 (oppb >> 3 & 0x1f)];
  421.         /* If they're not the same, dither between them. */
  422.         /* Dither constant is measured by where the true
  423.         color lies between the two nearest approximations.
  424.         Since the most nearly opposite color is not necessarily
  425.         on the line from the nearest through the true color,
  426.         some triangulation error can be introduced.  In the worst
  427.         case the r-nr distance can actually be less than the nr-oppr
  428.         distance. */
  429.         if (oppnrcm != nearestcm)
  430.         {
  431.         register int dither_const; 
  432.         int oppcmr = cmap[oppnrcm * 3];
  433.         int oppcmg = cmap[oppnrcm * 3 + 1];
  434.         int oppcmb = cmap[oppnrcm * 3 + 2];
  435.         dither_const = min(63, (64 * DIST(r, g, b, nr, ng, nb)) /
  436.                       DIST(nr, ng, nb, oppcmr, oppcmg, oppcmb));
  437.         if (!(dither_const >= 0 && dither_const < 64))
  438.         {
  439.             dbug2("Pix: %d, Dither constant: %d\n", pix, dither_const);
  440.             fprintf(stderr,"Near: %d %d %d\tTrue: %d %d %d\tOpp: %d %d %d / %d %d %
  441.             nr, ng, nb, r, g, b, oppr, oppg, oppb, oppcmr, oppcmg, oppcmb)
  442.         }
  443.         ASSERT (dither_const >= 0 && dither_const < 64);
  444.         if (Dither[(pix/cols)%DITH][(pix%cols)%DITH] < dither_const)
  445.             quant_image[pix] = oppnrcm;
  446.         else
  447.             quant_image[pix] = nearestcm;
  448.         } else
  449.         quant_image[pix] = nearestcm;
  450.     } else {
  451.         quant_image[pix] = nearestcm;
  452.     }
  453.     }
  454.     dbug0("Done.\n");
  455.     /* Cmap is still ok here */
  456. }
  457.  
  458. /* Load up the frame buffer color map - device specific */
  459. LoadMap(cmap, cmaplen)
  460. cmap_t *cmap;
  461. int cmaplen;
  462. {
  463.     register int i;
  464.     dbug1("Loading %d color map entries...", cmaplen);
  465. #ifdef GRAPHICS
  466.     for (i = 0; i < cmaplen*3; i+=3)
  467.     {
  468.     register unsigned int r = cmap[i], g = cmap[i+1], b = cmap[i+2];
  469.     mgicm(i/3, (r << 16) | (g << 8) | b);
  470.     dbug2("Color %d is 0x%x.\n", i/3, (r << 16) | (g << 8) | b);
  471.     }
  472. #endif
  473.     dbug0("done.\n");
  474. }
  475.  
  476. /* Draw the remapped image - device specific */
  477. Draw(image, rows, cols)
  478. out_pix_t *image;
  479. int rows, cols;
  480. {
  481. #ifdef GRAPHICS
  482.     mgiimage(image, 0, 0, cols-1, rows-1);
  483. #endif
  484. }
  485. SHAR_EOF
  486. fi
  487. if test -f 'split.c'
  488. then
  489.     echo shar: "will not over-write existing file 'split.c'"
  490. else
  491. cat << \SHAR_EOF > 'split.c'
  492. /*----------------------------------------------------------------------
  493.  * split up color space for quantize program
  494.  *----------------------------------------------------------------------
  495.  */
  496.  
  497. #include <stdio.h>
  498. #include <std.h>
  499. #include "quantize.h"
  500.  
  501. /* Pick your choice of algorithm: */
  502. /* MIDPOINT_CUT or POPULARITY */
  503. /* (Choose in makefile) */
  504.  
  505. #ifdef MIDPOINT_CUT
  506.  
  507. typedef unsigned int color_index_t;
  508.  
  509. /* My method here is to allocate a clist array big enough for all the
  510.  * possible colors, and then recursively split it.  Since the recursion
  511.  * is depth-first and we're building a LIST, not a TREE of boxes, we ca
  512.  * scramble the colors within a parent box because that parent won't be
  513.  * any more.  We end up with a list of boxes into which all the existin
  514.  * colors are partitioned, as described in Heckbert's paper.
  515.  */
  516. struct Box {
  517.     struct Box *next;
  518.     int rmin, rmax, gmin, gmax, bmin, bmax;
  519.     color_index_t *clist;    /* an array of color indices */
  520.     int clist_len;        /* how many cindexes in the list */
  521. };
  522.  
  523. static int nboxes, maxboxes;    /* number of boxes so far and max needed */
  524.  
  525. MakeCmap(hist, cmap, cmaplen)
  526. unsigned int *hist;
  527. cmap_t *cmap;
  528. int cmaplen;
  529. {
  530.     register int i;
  531.     struct Box *box;
  532.     nboxes = 1;
  533.     maxboxes = cmaplen;
  534.     NEW(box, struct Box, 1);
  535.     dbug0("Making color boxes.\n");
  536.     MakeTopBox(box, hist);    /* Make the first (top-level) box */
  537.     while (nboxes < maxboxes)
  538.         SplitBoxes(box);
  539.     FindColors(box, cmap, cmaplen);
  540. }
  541.  
  542. MakeTopBox(box, hist)        /* Make the top box from the histogram */
  543. struct Box *box;
  544. color_index_t *hist;
  545. {
  546.     register int i;
  547.     ASSERT(box != NULL);
  548.     ASSERT(hist != NULL);
  549.     
  550.     /* Allocate enough color cells to hold all the colors */
  551.     NEW(box->clist, color_index_t, HIST_SIZE);
  552.     box->clist_len = 0;
  553.     box->next = NULL;
  554.     /* For each possible rgb color */
  555.     for (i = 0; i < HIST_SIZE; i++)
  556.     {
  557.     if (hist[i] > 0)    /* if any pixels are that color, */
  558.         box->clist[box->clist_len++] = i;
  559.     }
  560.     Compact(box);
  561.     fprintf(stderr,"Top box at 0x%x: %d colors, R (%d -> %d)  G (%d -> %d)  B (%d -> %d)\n",
  562.         box, box->clist_len,
  563.         box->rmin, box->rmax, box->gmin, box->gmax, box->bmin, box->bmax);
  564. }
  565.  
  566. /* Shrink the box limits until it just encloses the points within it. */
  567. Compact(box)
  568. struct Box *box;
  569. {
  570.     register int i;
  571.     ASSERT(box != NULL);
  572.     box->rmin = 256;
  573.     box->gmin = 256;
  574.     box->bmin = 256;    /* Greater than any value */
  575.     box->rmax = -1;
  576.     box->gmax = -1;
  577.     box->bmax = -1;    /* Less than any rgb value */
  578.     for (i = box->clist_len-1; i >= 0; i--)
  579.     {
  580.     int r = (box->clist[i] >> 10) & 0x1f; /* colors associated with 'i' */
  581.     int g = (box->clist[i] >> 5) & 0x1f;
  582.     int b = (box->clist[i]) & 0x1f;
  583.     /* Expand the box until it just fits the colors tightly */
  584.     if (r < box->rmin) box->rmin = r;
  585.     if (g < box->gmin) box->gmin = g;
  586.     if (b < box->bmin) box->bmin = b;
  587.     if (r > box->rmax) box->rmax = r;
  588.     if (g > box->gmax) box->gmax = g;
  589.     if (b > box->bmax) box->bmax = b;
  590.     }
  591. }
  592.  
  593. /* Split all the boxes once each, stop if we've got enough */
  594. SplitBoxes(box)
  595. struct Box *box;
  596. {
  597.     struct Box *nextbox;
  598.     dbug1("Splitting all the boxes with %d done so far.\n", nboxes);
  599.     while (box != NULL && nboxes < maxboxes)
  600.     {
  601.     nextbox = box->next;
  602.     SplitBox(box);        /* SplitBox changes box->next! */
  603.     box = nextbox;
  604.     }
  605. }
  606.  
  607. static char longdim;        /* char representing the longest dimension */
  608.  
  609. int
  610. CompareColors(c1, c2)
  611. color_index_t *c1, *c2;
  612. {
  613.     switch (longdim)
  614.     {
  615.       case 'r':
  616.     return ((*c1 >> 10 & 0x1f) - (*c2 >> 10 & 0x1f));
  617.       case 'g':
  618.     return ((*c1 >> 5 & 0x1f) - (*c2 >> 5 & 0x1f));
  619.       case 'b':
  620.     return ((*c1 & 0x1f) - (*c2 & 0x1f));
  621.       default:
  622.     fprintf(stderr,"Illegal longdim 0x%x in CompareColors.\n", longdim);
  623.     exit(1);
  624.     }
  625. }
  626.  
  627. /* Divide the box across its longest dimension at its median point
  628.  * This is NOT recursive - we split all the boxes once, then if there's
  629.  * more colors to be assigned we split all of them again, and so on. */
  630. SplitBox(box)
  631. struct Box *box;
  632. {
  633.     struct Box *newbox;        /* the 'other half' of the split box */
  634.     ASSERT(box != NULL);
  635.     if (nboxes > maxboxes) return;
  636.     if (box->clist_len <= 1) return;
  637.     nboxes++;
  638.     if (box->rmax - box->rmin >= box->gmax - box->gmin &&
  639.     box->rmax - box->rmin >= box->bmax - box->bmin)
  640.     longdim = 'r';
  641.     else
  642.     if (box->gmax - box->gmin >= box->rmax - box->rmin &&
  643.     box->gmax - box->gmin >= box->bmax - box->bmin)
  644.     longdim = 'g';
  645.     else
  646.     if (box->bmax - box->bmin >= box->gmax - box->gmin &&
  647.     box->bmax - box->bmin >= box->rmax - box->rmin)
  648.     longdim = 'b';
  649.     else
  650.     {
  651.     fprintf(stderr,"No longest dimension in SplitBox??!\n");
  652.         fprintf(stderr,"This box: R (%d -> %d)  G (%d -> %d)  B (%d -> %d)\n",
  653.         box->rmin, box->rmax, box->gmin, box->gmax, box->bmin, box-
  654.     exit(1);
  655.     }
  656.     qsort(box->clist, box->clist_len, sizeof(color_index_t), CompareColors);
  657.     /* Set up the new box and put it after the old box */
  658.     NEW(newbox, struct Box, 1);
  659.     newbox->next = box->next;
  660.     box->next = newbox;
  661.     /* If len is odd, put the extra one in the old box */
  662.     newbox->clist_len = box->clist_len / 2;
  663.     box->clist_len -= newbox->clist_len;
  664.     newbox->clist = box->clist + box->clist_len;
  665.     Compact(box);
  666.     Compact(newbox);
  667. }
  668.  
  669. /* Find the representative color in each box and put it into cmap. */
  670. FindColors(box, cmap, cmaplen)    /* here box is a list of boxes */
  671. struct Box *box;
  672. cmap_t *cmap;
  673. int cmaplen;
  674. {
  675.     register int r, g, b;
  676.     register int i;
  677.     int cml = 0;
  678.     while (box != NULL)
  679.     {
  680.     r = 0, g = 0, b = 0;
  681.     for (i=0; i < box->clist_len; i++)
  682.     {
  683.         r += (box->clist[i] >> 10 & 0x1f) << 3;
  684.         g += (box->clist[i] >> 5 & 0x1f) << 3;
  685.         b += (box->clist[i] & 0x1f) << 3;
  686.     }
  687.     ASSERT(r/box->clist_len < 256);
  688.     ASSERT(g/box->clist_len < 256);
  689.     ASSERT(b/box->clist_len < 256);
  690.     cmap[cml*3] = r / box->clist_len;
  691.     cmap[cml*3+1] = g / box->clist_len;
  692.     cmap[cml*3+2] = b / box->clist_len;
  693.     box = box->next;
  694.     cml++;
  695.     }
  696.     ASSERT(cmaplen == cml);    /* since there should be 'cmaplen' boxes */
  697.     for (i = 0; i < cmaplen*3; i+=3)
  698.     {
  699.     printf("Cmap[%4d] = %3d, %3d, %3d\n",
  700.        i/3, cmap[i], cmap[i+1], cmap[i+2]);
  701.     }
  702. }
  703.  
  704. #endif
  705.  
  706. #ifdef POPULARITY
  707. static unsigned int *HiSt;
  708.  
  709. CompareHist(v1, v2)
  710. int *v1, *v2;
  711. {
  712.     return (HiSt[*v2] - HiSt[*v1]);
  713. }
  714.  
  715. MakeCmap(hist, cmap, cmaplen)
  716. unsigned int *hist;
  717. cmap_t *cmap;
  718. int cmaplen;
  719. {
  720.     register int i;
  721.     unsigned int *hptrs;
  722.     /* This is the popularity algorithm (Boyle & Lippman, Archmac, 1978) */
  723.     /* First we sort the histogram: set up a bunch of pointers and 
  724.     dbug1("Getting %d ints for hptrs...", HIST_SIZE);
  725. /*  NEW(hptrs, unsigned int, HIST_SIZE);*/
  726.     hptrs = (unsigned int *)quant_image;    /* use existing storage temporarily */
  727.     dbug0("initting it...");
  728.     for (i = 0; i < HIST_SIZE; i++) hptrs[i] = i;
  729.     HiSt = hist;
  730.     dbug0("Sorting...");
  731.     qsort (hptrs, HIST_SIZE, sizeof(unsigned int), CompareHist);
  732.     dbug0("making color map...");
  733.     for (i = 0; i < cmaplen; i++)
  734.     {
  735.     cmap[3*i] = hptrs[i] >> 10 << 3;
  736.     cmap[3*i+1] = (hptrs[i] >> 5 & 0x1f) << 3;
  737.     cmap[3*i+2] = (hptrs[i] & 0x1f) << 3;
  738.     }
  739.     dbug0("done.\n");
  740.     dbug1("Max pixels in one bucket: %d.\n", hist[hptrs[0]]);
  741. }
  742. #endif POPULARITY
  743.  
  744. SHAR_EOF
  745. fi
  746. if test -f 'quantize.h'
  747. then
  748.     echo shar: "will not over-write existing file 'quantize.h'"
  749. else
  750. cat << \SHAR_EOF > 'quantize.h'
  751. /*----------------------------------------------------------------------
  752.  * header file for color-quantizer
  753.  * by Gary Oberbrunner  1/5/88
  754.  *----------------------------------------------------------------------
  755.  */
  756.  
  757. #define HIST_QUANT_BITS 5    /* Histogram clump quantization */
  758. #define HIST_SHIFT (8 - HIST_QUANT_BITS)
  759. #define HIST_SIZE (1 << (HIST_QUANT_BITS * 3))
  760.  
  761. #define RGB_BYTES 4        /* Bytes per rgb pixel (4 for alpha channel) */
  762.  
  763. /* We need the next definition to reference a variable-size array */
  764. /* NCOLS is the number of columns in the array */
  765. /* Each pixel is stored as an rgb byte triplet */
  766.  
  767. /* This is the address of the pixel triplet: */
  768. #define RCPIX(array, row, col)    (array + RGB_BYTES * (row * NCOLS + col))
  769. #define PIXPIX(array, pixel)    (array + RGB_BYTES * (pixel))
  770. /* These are the pixel components: */
  771. #define RED(array, row, col)    (*(RCPIX(array, row, col)))
  772. #define GREEN(array, row, col)    (*(RCPIX(array, row, col) + 1))
  773. #define BLUE(array, row, col)    (*(RCPIX(array, row, col) + 2))
  774. #define PIXRED(array, pix)    (*(PIXPIX(array, pix)))
  775. #define PIXGREEN(array, pix)    (*(PIXPIX(array, pix) + 1))
  776. #define PIXBLUE(array, pix)    (*(PIXPIX(array, pix) + 2))
  777.  
  778. /* Here is the rgb-space distance metric: */
  779. #define DIST(r1,g1,b1,r2,g2,b2) \
  780.        (int) \
  781.        (3 * ((r1)-(r2)) * ((r1)-(r2)) + \
  782.     4 * ((g1)-(g2)) * ((g1)-(g2)) + \
  783.     2 * ((b1)-(b2)) * ((b1)-(b2)))
  784. #define BIG_DISTANCE 1000000    /* bigger than any real distance */
  785.  
  786. typedef unsigned char cmap_t;
  787. typedef unsigned char rgb_pix_t;
  788. typedef unsigned char out_pix_t;
  789.  
  790. extern out_pix_t  *quant_image;    /* the quantized image - cmap indices */
  791. extern int debug;
  792. SHAR_EOF
  793. fi
  794. if test -f 'std.h'
  795. then
  796.     echo shar: "will not over-write existing file 'std.h'"
  797. else
  798. cat << \SHAR_EOF > 'std.h'
  799. /* Standard include file with lots of generally useful stuff */
  800. /* by Gary Oberbrunner */
  801.  
  802. #define min(x,y) ((x) < (y) ? (x) : (y))
  803. #define max(x,y) ((x) > (y) ? (x) : (y))
  804. #define clamp(x, mn, mx) (((x) <= (mn)) ? (mn) : (((x) >= (mx)) ? (mx) : (x)))
  805.  
  806. #define SECURITY
  807. #ifdef SECURITY
  808. #define ASSERT(exp) { if (!(exp)) { fprintf(stderr,\
  809.     "Assertion error: %s, line %d. Assertion exp failed.\n",\
  810.     __FILE__, __LINE__);\
  811.     exit(69); } }
  812. #else
  813. #define ASSERT(exp)
  814. #endif
  815.  
  816. #define dbug0(str) \
  817.     if (debug) fprintf(stderr, str);
  818. #define dbug1(str, arg1) \
  819.     if (debug) fprintf(stderr, str, arg1);
  820. #define dbug2(str, arg1, arg2) \
  821.     if (debug) fprintf(stderr, str, arg1, arg2);
  822. #define dbug3(str, arg1, arg2, arg3) \
  823.     if (debug) fprintf(stderr, str, arg1, arg2, arg3);
  824. #define dbug4(str, arg1, arg2, arg3, arg4) \
  825.     if (debug) fprintf(stderr, str, arg1, arg2, arg3, arg4);
  826.  
  827. extern char *malloc();
  828. /*
  829.  * NEW is a macro to malloc 'n' new variables of type 'type'.
  830.  */
  831. #define NEW(var, type, num) \
  832.     if ((var = (type *) malloc((num) * sizeof(type)))==NULL) \
  833.     {    fprintf(stderr,\
  834.         "ERROR: Out of memory.\nLast request:\n");\
  835.     fprintf(stderr,\
  836.         "\tNumber: num\n\tType  : type\n\tName  : var.\n");\
  837.     fprintf(stderr,"File %s, line %d\n", __FILE__, __LINE__);\
  838.     fprintf(stderr,\
  839.     "\tTotal requested: num=%d x sizeof(type)=%d\n\t =%d bytes.\n",\
  840.         num, sizeof(type), (num) * sizeof(type));\
  841.     exit(99); }
  842. SHAR_EOF
  843. fi
  844. exit 0
  845. #    End of shell archive
  846. -- 
  847. Remember,               -Truth is not beauty;
  848. Information is not knowledge; /    Beauty is not love;      Gary Oberbrunner
  849. Knowledge is not wisdom;     /    Love is not music;      ...!masscomp!garyo
  850. Wisdom is not truth;    ----/    Music is the best. - FZ
  851. -----------------------------------------------------------------------------
  852. From: awpaeth@watcgl.waterloo.edu (Alan W. Paeth)
  853. Subject: Re: Advanced Dither Needed
  854. Date: 8 Jan 88 22:41:45 GMT
  855.  
  856. A comment on dither - the "quick and dirty" approach to uniform division of
  857. the color space (so that R, G and B can be treated separably) very 
  858. slices up 8 bit pixels as 3 bits for red and green, 2 bits for blue (where
  859. the eye is least sensitive).
  860.  
  861. This is an unnecessary oversimplification that leaves blue with only two
  862. mid-range intensities. It suggests itself because we tend to think in 
  863. of "bits" per color at the hardware level, not "discrete number of inten
  864. steps". This approach is further suggested because the pixel channels 
  865. R G and B can be written separably through the proper setting of the wri
  866. mask, but in practice, the pixel is normally written as a single byte 
  867. not as three consecutive "inking" passes. One useful property of this ap
  868. is that there are a number of "greys" in the color space, as the chann
  869. precisions are all multiples of each other.
  870.  
  871. A better approach is to form a color table containing a Cartesian product of
  872. intensites on a range that is not necessarily a power of two. For 
  873. we take red = [0..5], green = [0..6], blue = [0..4] (suitably normalized to
  874. represent intensities on the range [0.0..1.0] then there are 6*7*5 
  875. table entries, and noticibly better representation in blue.
  876.  
  877. Our Graphics Lab takes this approach for writing on our VAX GPX's under X.
  878. We must use a uniform space, because many simultaneous image (window
  879. be on the display, and each would have a different histogram (if Heckbert'
  880. approach were used), potentially overflowing the color table. In our
  881. X-write IM tool takes color images of arbitrary precision and Floyd-Steinb
  882. converts them into this space of 6,7,5 steps, and thus all color ima
  883. identical entries within the color table (which is a consequence of how X
  884. manages requests for color table entries). B/W images typically alloc
  885. level grey ramp which can co-exist within the remaining space, leaving
  886. 14 (256=210+32+14) entries for use by other applications, should they re
  887. colors other than saturated yellow, blue, white, etc.
  888.  
  889. On GPX's with smaller color tables, the same approach is used. Once again, ~3/4
  890. of the hardware color table is allocated, giving 12 free color 
  891. the program divides as 2*3*2 in RGB. Other useful sets have been 6*6*6, which
  892. gives the largest uniform allocation in 256 entries, thus allowin
  893. greys, and 6*7*6, which uses 252 out of 256 colors and is much better overall
  894. than the 8*8*4 set which is used implicitly when one allocates us
  895. "three bits red and green, two for blue" scheme.
  896.  
  897.     /Alan Paeth
  898.     Computer Graphics Lab
  899.     University of Waterloo
  900. -----------------------------------------------------------------------------
  901. From: lou@aramis.rutgers.edu (Lou Steinberg)
  902. Subject: Re: Advanced Dither Needed
  903. Date: 11 Jan 88 04:46:33 GMT
  904.  
  905. In article <2741@masscomp.UUCP> garyo@masscomp.UUCP (Gary Oberbrunner) writes:
  906.  
  907. > [for] Floyd & Steinberg's dithering algorithm [...]
  908. > On some images it can help to add some random noise; for instance if you
  909. > have a smooth (flat-shaded) surface that's just below one of the c
  910. > the map, the error will build up slowly until you get a line of brighter
  911. > color somewhere in the surface.  A bit o' noise gets rid of that l
  912.  
  913. The other way to get rid of it is to alternate scanning left to right
  914. and right to left, i.e. even scan lines one way, odd scan lines the
  915. other.  The line of brighter color essentially comes because pixels on
  916. that line do not communicate - there is no path for the error value at
  917. one to affect the value chosen for the other.  By scanning in a zig
  918. zag fashion, you ensure that each pixel communicates with every other
  919. pixel scanned after it.
  920. -- 
  921.                     Lou Steinberg
  922.  
  923. uucp:   {pretty much any major site}!rutgers!aramis.rutgers.edu!lou 
  924. arpa:   lou@aramis.rutgers.edu
  925. -----------------------------------------------------------------------------
  926. From: lou@aramis.rutgers.edu (Lou Steinberg)
  927. Subject: Re: Advanced Dither Needed
  928. Date: 11 Jan 88 05:02:13 GMT
  929.  
  930. In article <2741@masscomp.UUCP> garyo@masscomp.UUCP (Gary Oberbrunner) writes:
  931. > But the wayd [Floyd/Steinberg] dither works is like this (it's
  932.   [...]
  933. > Add 3/8 of that error to the pixel to the right, and 3/8 to the pixel
  934. > below the current pixel.  The remaining 1/4 goes to the pixel diagona
  935. > below.
  936.  
  937. I should point out that the actual fractions we used were, assuming
  938. you are at X, moving left to right:
  939.  
  940.         X     7/16
  941.          3/16   5/16  1/16    
  942.  
  943. Note that the error goes to four neighbors, not three.  I think this
  944. will probably do better (at least for black and white) than the
  945. 3/8-3/8-1/4 distribution, at the cost of greater processing.  I have
  946. seen the 3/8-3/8-1/4 distribution described as "our" algorithm before,
  947. but I have no idea who the credit really belongs to.
  948.  
  949. Also, I should add that if you do zig-zag scanning (see my immediately
  950. previous message), it is sufficient (but not quite as good) to send
  951. half the error one pixel ahead (e.g. to the right on lines you scan
  952. left to right), and half one pixel straight down.  Again, this is for
  953. black and white;  I've not tried it with color.
  954. -- 
  955.                     Lou Steinberg
  956.  
  957. uucp:   {pretty much any major site}!rutgers!aramis.rutgers.edu!lou 
  958. arpa:   lou@aramis.rutgers.edu
  959.  
  960. Note: In a seperate message, another person confirmed that this method
  961.       works well with color too.
  962.  
  963.